{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spatial derivatives of $\\lambda$ for oblate ellipsoids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Semi-axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = 20.\n",
    "b = 1.3*a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random points $(x, y, z)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 1000\n",
    "R = 1.1\n",
    "theta, phi = np.meshgrid(np.linspace(0., np.pi, N), np.linspace(0., 2.*np.pi, N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = a*R*np.cos(theta)\n",
    "y = b*R*np.sin(theta)*np.cos(phi)\n",
    "z = b*R*np.sin(theta)*np.sin(phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coefficients of the quadratic function for calculating $\\lambda$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p2 = 1.\n",
    "p1 = a**2 + b**2 - (x**2) - (y**2) - (z**2)\n",
    "p0 = (a*b)**2 - (b*x)**2 - (a*y)**2 - (a*z)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameter $\\lambda$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "delta = p1**2 - 4.*p2*p0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lamb = (-p1 + np.sqrt(delta))/(2.*p2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spatial derivatives calculated here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dlamb_dx = (2.*x/(a**2 + lamb))/((x/(a**2 + lamb))**2 + (y/(b**2 + lamb))**2 + (z/(b**2 + lamb))**2)\n",
    "dlamb_dy = (2.*y/(b**2 + lamb))/((x/(a**2 + lamb))**2 + (y/(b**2 + lamb))**2 + (z/(b**2 + lamb))**2)\n",
    "dlamb_dz = (2.*z/(b**2 + lamb))/((x/(a**2 + lamb))**2 + (y/(b**2 + lamb))**2 + (z/(b**2 + lamb))**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spatial derivatives presented by Emerson et al. (1985)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r2 = x**2 + y**2 + z**2\n",
    "delta = np.sqrt(r2*r2 + (a**2 - b**2)**2 - 2.*(a**2 - b**2)*(x**2 - y**2 - z**2))\n",
    "\n",
    "dlamb_dx_Emerson85 = x*(1. + (r2 - a**2 + b**2)/delta)\n",
    "dlamb_dy_Emerson85 = y*(1. + (r2 + a**2 - b**2)/delta)\n",
    "dlamb_dz_Emerson85 = z*(1. + (r2 + a**2 - b**2)/delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison between the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(dlamb_dx, dlamb_dx_Emerson85)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(dlamb_dy, dlamb_dy_Emerson85)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(dlamb_dz, dlamb_dz_Emerson85)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
